COVID-19 US

Prepare county level data

Read and format personality data

df_us_pers <- read_csv('timeseries_usa_county_march1_april_09.csv')

df_us_pers <- df_us_pers %>% select(countyfips, open, sci, extra, agree, stabil) %>% 
  mutate(stabil = 6-stabil) %>%
  dplyr::rename(county_fips = countyfips,
         pers_o = open, 
         pers_c = sci,
         pers_e = extra,
         pers_a = agree,
         pers_n = stabil) %>% 
  distinct() %>%
  mutate(county_fips = as.character(county_fips))

Read and format prevalence data

df_us_prev <- read_csv('USA_timeseries_prep_2005.csv')

df_us_prev <- df_us_prev %>% 
  select(fips, date, rate) %>% 
  mutate(date = as.Date(date, "%d%b%Y")) %>% 
  rename(county_fips = fips, 
         rate_day = rate) %>%
  mutate(county_fips = as.character(county_fips))

df_us_prev

Read and format county level controls

df_us_ctrl <- read.csv('controls_US.csv')

df_us_ctrl <- df_us_ctrl %>% select(-county_name) %>% 
  rename(county_fips = county) %>%
  mutate(county_fips = as.character(county_fips))

df_us_ctrl

Read and format social distancing data FB

fb_files <- list.files('../FB Data/US individual files/Mobility/',
                       '*.csv', full.names = T)

df_us_socdist <- fb_files %>% 
  map(read_csv) %>% bind_rows()

df_us_socdist <- df_us_socdist %>%
  select(-age_bracket, -gender, -baseline_name, -baseline_type, -polygon_name) %>%
  rename(date = ds,
         county_fips = polygon_id,
         socdist_tiles = all_day_bing_tiles_visited_relative_change,
         socdist_single_tile = all_day_ratio_single_tile_users) %>%
  mutate(county_fips = as.character(county_fips))

df_us_socdist

df_us_socdist %>% select(county_fips) %>% distinct() %>% nrow()

Merge data

# create sequence of dates
date_sequence <- seq.Date(min(df_us_socdist$date),
                          as.Date('2020-04-28'), 1)
                     
# create data frame with time sequence
df_dates = data.frame(date_sequence, 1:length(date_sequence)) 
names(df_dates) <- c('date', 'time')

# join data frames 
df_us_socdist <- df_us_socdist %>%
  plyr::join(df_us_ctrl, by='county_fips') %>% 
  plyr::join(df_us_pers, by='county_fips') %>%
  inner_join(df_dates, by='date') %>% 
  arrange(county_fips, date)

fips_complete <- df_us_socdist %>% 
  group_by(county_fips) %>% 
  summarize(n = n()) %>% 
  filter(! n<max(n)) %>% .$county_fips

df_us_socdist <- df_us_socdist %>%
  filter(county_fips %in% fips_complete)

df_us_socdist

Control for weekend effect

easter <- seq.Date(as.Date('2020-04-10'), as.Date('2020-04-13'), 1)

df_us_loess <- df_us_socdist %>% 
  mutate(weekday = format(date, '%u')) %>% 
  filter(!weekday %in% c('6','7') | date %in% easter) %>% 
  split(.$county_fips) %>%
  map(~ loess(socdist_single_tile ~ time, data = .)) %>%
  map(predict, 1:max(df_us_socdist$time)) %>% 
  bind_rows() %>% 
  gather(key = 'county_fips', value = 'loess') %>% 
  group_by(county_fips) %>% 
  mutate(time = row_number())

df_us_loess_2 <- df_us_socdist %>% 
  mutate(weekday = format(date, '%u')) %>% 
  filter(!weekday %in% c('6','7') | date %in% easter) %>% 
  split(.$county_fips) %>%
  map(~ loess(socdist_tiles ~ time, data = .)) %>%
  map(predict, 1:max(df_us_socdist$time)) %>% 
  bind_rows() %>% 
  gather(key = 'county_fips', value = 'loess') %>% 
  rename(loess_2 = loess) %>%
  group_by(county_fips) %>% 
  mutate(time = row_number())

df_us_socdist <- df_us_socdist %>% 
  merge(df_us_loess, by=c('county_fips', 'time')) %>% 
  merge(df_us_loess_2, by=c('county_fips', 'time')) %>% 
  mutate(weekday = format(date, '%u')) %>% 
  mutate(socdist_single_tile_clean = ifelse(weekday %in% c('6','7') | date %in% easter, 
                                            loess, socdist_single_tile),
         socdist_tiles_clean = ifelse(weekday %in% c('6','7') | date %in% easter, 
                                            loess_2, socdist_tiles)) %>%
  arrange(county_fips, time) %>% 
  select(-weekday)

df_us_socdist <- df_us_socdist %>% drop_na() %>% mutate(time = time-1)

Plot prevalence over time

df_us_prev %>% sample_n(20000) %>%
  ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=county_fips, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall prevalence over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_us_prev %>% 
  mutate(prev_tail = cut(.[[i]], 
                         breaks = c(-Inf, quantile(.[[i]], 0.05, na.rm=T), 
                                    quantile(.[[i]], 0.95, na.rm=T), Inf),
                         labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(prev_tail != 'center') %>%
  ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=county_fips, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~prev_tail) + 
  theme(legend.position="none") + 
  ggtitle(i)

print(gg)
}

Plot social distancing single tile visited

df_us_socdist %>% sample_n(10000) %>%
  ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=county_fips, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall social distancing (single tile) over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_us_socdist %>% 
  mutate(dist_tail = cut(.[[i]], 
                         breaks = c(-Inf, quantile(.[[i]], 0.05, na.rm = T), 
                                    quantile(.[[i]], 0.95, na.rm = T), Inf),
                         labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(dist_tail != 'center') %>%
  ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=county_fips, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~dist_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

df_us_socdist %>% sample_n(10000) %>%
  ggplot(aes(x=time, y=socdist_tiles_clean)) + 
  geom_point(aes(col=county_fips, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall social distancing (single tile) over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_us_socdist %>% 
  mutate(dist_tail = cut(.[[i]], 
                         breaks = c(-Inf, quantile(.[[i]], 0.05, na.rm = T), 
                                    quantile(.[[i]], 0.95, na.rm = T), Inf),
                         labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(dist_tail != 'center') %>%
  ggplot(aes(x=time, y=socdist_tiles_clean)) + 
  geom_point(aes(col=county_fips, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~dist_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

Example plots

df_us_socdist <- df_us_socdist %>% 
  mutate(socdist_single_tile = socdist_single_tile_clean,
         socdist_tiles = socdist_tiles_clean) %>% 
  select(-loess, -loess_2, -socdist_single_tile_clean, -socdist_tiles_clean)

df_us_prev %>% select(-time, -date) %>% 
  group_by(county_fips) %>%
  summarize_if(is.numeric, mean) %>% 
  select(-county_fips) %>%
  cor(use='pairwise.complete.obs') %>% 
  round(3) %>% as.data.frame()

df_us_socdist %>% select(-time, -date) %>% 
  group_by(county_fips) %>%
  summarize_if(is.numeric, mean) %>% 
  select(-county_fips) %>%
  cor(use='pairwise.complete.obs') %>% 
  round(3) %>% as.data.frame()

NA

lvl2_scaled <- df_us_prev %>% 
  select(-time, -date, -rate_day) %>% 
  distinct() %>% 
  mutate_at(vars(-county_fips), scale)

lvl1_scaled <- df_us_prev %>% select(county_fips, time, rate_day)

df_us_prev_scaled <- plyr::join(lvl1_scaled, lvl2_scaled, by = 'county_fips') 

df_us_socdist <- df_us_socdist %>% 
  mutate(socdist_single_tile = socdist_single_tile_clean,
         socdist_tiles = socdist_tiles_clean) %>% 
  select(-loess, -loess_2, -socdist_single_tile_clean, -socdist_tiles_clean)

Correlations

lvl2_scaled <- df_us_socdist %>% 
  select(-time, -date, -socdist_tiles, -socdist_single_tile) %>% 
  distinct() %>% 
  mutate_at(vars(-county_fips), scale)

lvl1_scaled <- df_us_socdist %>% 
  select(county_fips, time, socdist_single_tile, socdist_tiles)

df_us_socdist_scaled <- plyr::join(lvl1_scaled, lvl2_scaled, by = 'county_fips') 

Rescale Data

# get onset day
df_us_onset_prev <- df_us_prev_scaled %>% 
  group_by(county_fips) %>% 
  mutate(rate_cs = cumsum(rate_day)) %>% 
  filter(rate_cs > 0) %>%
  summarize(onset_prev = min(time)) %>%
  mutate(county_fips = as.character(county_fips))
  
# merge with county data
df_us_onset_prev <- df_us_prev_scaled %>% 
  select(-time, -rate_day) %>%
  distinct() %>% 
  mutate(county_fips = as.character(county_fips)) %>%
  left_join(df_us_onset_prev, by = 'county_fips')

# handle censored data
df_us_onset_prev <- df_us_onset_prev %>% 
  mutate(event = ifelse(is.na(onset_prev), 0, 1)) %>% 
  mutate(onset_prev = replace_na(onset_prev, as.numeric(diff(range(df_us_prev$date)))+1))

lvl2_scaled <- df_us_socdist %>% 
  select(-time, -date, -socdist_tiles, -socdist_single_tile) %>% 
  distinct() %>% 
  mutate_at(vars(-county_fips), scale)

lvl1_scaled <- df_us_socdist %>% 
  select(county_fips, time, socdist_single_tile, socdist_tiles)

df_us_socdist_scaled <- plyr::join(lvl1_scaled, lvl2_scaled, by = 'county_fips') 

Predict Prevalence

Extract first day of covid outbreak

# get onset day
df_us_onset_prev <- df_us_prev_scaled %>% 
  group_by(county_fips) %>% 
  mutate(rate_cs = cumsum(rate_day)) %>% 
  filter(rate_cs > 0) %>%
  summarize(onset_prev = min(time)) %>%
  mutate(county_fips = as.character(county_fips))
  
# merge with county data
df_us_onset_prev <- df_us_prev_scaled %>% 
  select(-time, -rate_day) %>%
  distinct() %>% 
  mutate(county_fips = as.character(county_fips)) %>%
  left_join(df_us_onset_prev, by = 'county_fips')

# handle censored data
df_us_onset_prev <- df_us_onset_prev %>% 
  mutate(event = ifelse(is.na(onset_prev), 0, 1)) %>% 
  mutate(onset_prev = replace_na(onset_prev, as.numeric(diff(range(df_us_prev$date)))+1))

Extract slopes

# cut time series before onset
df_us_prev_scaled <- df_us_prev_scaled %>% 
  group_by(county_fips) %>% 
  mutate(rate_cs = cumsum(rate_day)) %>% 
  filter(rate_cs > 0) %>%
  mutate(time = time-min(time)+1) %>%
  ungroup() %>%
  filter(time <= 30) %>%
  select(-rate_cs)

# drop counties with little data
df_us_prev_scaled <- df_us_prev_scaled %>%
  group_by(county_fips) %>%
  filter(n() == 30) %>%
  ungroup()

# log transform prevalence data 
df_us_prev_scaled <- df_us_prev_scaled %>% 
  mutate(rate_day = log(rate_day)) %>%
  filter(rate_day != -Inf)

# extract slope prevalence
df_us_slope_prev <- df_us_prev_scaled %>% split(.$county_fips) %>% 
  map(~ lm(rate_day ~ time, data = .)) %>%
  map(coef) %>% 
  map_dbl('time') %>% 
  as.data.frame() %>% 
  rownames_to_column('county_fips') %>% 
  rename(slope_prev = '.')
  
# merge with county data
df_us_slope_prev <- df_us_onset_prev %>% 
  inner_join(df_us_slope_prev, by = 'county_fips') %>%
  drop_na()

# handle censored data
df_us_slope_prev_cens <- df_us_slope_prev %>% 
  filter(slope_prev>=0) %>% 
  mutate(slope_prev = scale(slope_prev),
         onset_prev = scale(onset_prev))

# standardize slopes
df_us_slope_prev <- df_us_slope_prev %>% 
  mutate(slope_prev = scale(slope_prev),
         onset_prev = scale(onset_prev))

Explore distributions

# predict slopes from personality
lm_slope_prev <- lm(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_us_slope_prev)
lm_slope_prev %>% summary()

Call:
lm(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n, data = df_us_slope_prev)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.8951 -0.6401 -0.0422  0.5718  3.4926 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.01667    0.01996  -0.835    0.404    
pers_o       0.26702    0.02277  11.726  < 2e-16 ***
pers_c       0.01018    0.02759   0.369    0.712    
pers_e       0.09140    0.02290   3.992 6.76e-05 ***
pers_a       0.02662    0.02812   0.947    0.344    
pers_n      -0.11496    0.02584  -4.450 9.02e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9472 on 2253 degrees of freedom
Multiple R-squared:  0.1048,    Adjusted R-squared:  0.1028 
F-statistic: 52.75 on 5 and 2253 DF,  p-value: < 2.2e-16

lm_slope_prev %>% confint(level=0.9)

                    5 %        95 %
(Intercept) -0.04951426  0.01617791
pers_o       0.22954727  0.30449288
pers_c      -0.03522143  0.05558882
pers_e       0.05372148  0.12907279
pers_a      -0.01965013  0.07289248
pers_n      -0.15747983 -0.07244929

# predict slopes from personality with controls
lm_slope_prev_ctrl <- lm(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n + 
                           airport_distance + republican + medage + male + popdens + 
                           manufact + tourism + academics + medinc + physician_pc +
                           onset_prev,
                         data = df_us_slope_prev)
lm_slope_prev_ctrl %>% summary()

Call:
lm(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + airport_distance + republican + medage + male + 
    popdens + manufact + tourism + academics + medinc + physician_pc + 
    onset_prev, data = df_us_slope_prev)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.7005 -0.4837 -0.0430  0.4396  4.3316 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -0.068070   0.017031  -3.997 6.63e-05 ***
pers_o            0.037036   0.021795   1.699  0.08939 .  
pers_c            0.047875   0.022496   2.128  0.03343 *  
pers_e            0.022047   0.018639   1.183  0.23699    
pers_a            0.046081   0.023623   1.951  0.05122 .  
pers_n            0.031360   0.022704   1.381  0.16733    
airport_distance -0.135510   0.022505  -6.021 2.02e-09 ***
republican       -0.118138   0.020286  -5.824 6.59e-09 ***
medage           -0.133828   0.018098  -7.395 1.99e-13 ***
male             -0.049010   0.021397  -2.291  0.02208 *  
popdens           0.096092   0.016690   5.757 9.71e-09 ***
manufact         -0.006815   0.019549  -0.349  0.72739    
tourism           0.012977   0.021499   0.604  0.54617    
academics        -0.142904   0.030820  -4.637 3.74e-06 ***
medinc            0.183557   0.026380   6.958 4.51e-12 ***
physician_pc     -0.047652   0.017896  -2.663  0.00781 ** 
onset_prev       -0.428032   0.020786 -20.592  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7568 on 2242 degrees of freedom
Multiple R-squared:  0.4314,    Adjusted R-squared:  0.4273 
F-statistic: 106.3 on 16 and 2242 DF,  p-value: < 2.2e-16

lm_slope_prev_ctrl %>% confint(level=0.9)

                          5 %        95 %
(Intercept)      -0.096094777 -0.04004524
pers_o            0.001172681  0.07289998
pers_c            0.010857177  0.08489190
pers_e           -0.008623806  0.05271706
pers_a            0.007207823  0.08495406
pers_n           -0.005999178  0.06871970
airport_distance -0.172542458 -0.09847698
republican       -0.151519389 -0.08475717
medage           -0.163608583 -0.10404834
male             -0.084219727 -0.01380084
popdens           0.068628060  0.12355653
manufact         -0.038983249  0.02535238
tourism          -0.022400681  0.04835494
academics        -0.193619372 -0.09218857
medinc            0.140148361  0.22696623
physician_pc     -0.077099345 -0.01820385
onset_prev       -0.462237042 -0.39382737

Predict COVID onset with time-to-event regression

# predict onset from personality
cox_onset_prev <- coxph(Surv(onset_prev, event) ~ 
                          pers_o + pers_c + pers_e + pers_a + pers_n, 
                        data = df_us_onset_prev)
cox_onset_prev %>% summary()

# predict onset from personality with controls
cox_onset_prev_ctrl <- coxph(Surv(onset_prev, event) ~ 
                               pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_distance + republican + medage + male + popdens +
                               manufact + tourism + academics + medinc + physician_pc,
                             data = df_us_onset_prev)
cox_onset_prev_ctrl %>% summary()

Predict prevalence slopes with linear models

# predict slopes from personality
tobit_slope_prev <- censReg(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n, 
                            left = min(df_us_slope_prev_cens$slope_prev), right = Inf,
                         data = df_us_slope_prev_cens)
tobit_slope_prev %>% summary()

Call:
censReg(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n, left = min(df_us_slope_prev_cens$slope_prev), right = Inf, 
    data = df_us_slope_prev_cens)

Observations:
         Total  Left-censored     Uncensored Right-censored 
          2237              6           2231              0 

Coefficients:
            Estimate Std. error t value  Pr(> t)    
(Intercept) -0.01884    0.02006  -0.939 0.347772    
pers_o       0.27674    0.02296  12.054  < 2e-16 ***
pers_c       0.01671    0.02778   0.601 0.547523    
pers_e       0.10728    0.02319   4.625 3.74e-06 ***
pers_a       0.01551    0.02845   0.545 0.585543    
pers_n      -0.10392    0.02607  -3.986 6.72e-05 ***
logSigma    -0.05453    0.01498  -3.640 0.000273 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Newton-Raphson maximisation, 3 iterations
Return code 1: gradient close to zero
Log-likelihood: -3052.156 on 7 Df

tobit_slope_prev %>% confint(level=0.9)

                    5 %        95 %
(Intercept) -0.05183356  0.01416185
pers_o       0.23897737  0.31450381
pers_c      -0.02898477  0.06240249
pers_e       0.06912983  0.14542999
pers_a      -0.03128303  0.06231182
pers_n      -0.14680511 -0.06103538
logSigma    -0.07917392 -0.02988752

# predict slopes from personality with controls
tobit_slope_prev_ctrl <- censReg(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n + 
                           airport_distance + republican + medage + male + popdens + 
                           manufact + tourism + academics + medinc + physician_pc +
                             onset_prev,
                           left = min(df_us_slope_prev_cens$slope_prev), right = Inf,
                         data = df_us_slope_prev_cens)
tobit_slope_prev_ctrl %>% summary()

Call:
censReg(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + airport_distance + republican + medage + male + 
    popdens + manufact + tourism + academics + medinc + physician_pc + 
    onset_prev, left = min(df_us_slope_prev_cens$slope_prev), 
    right = Inf, data = df_us_slope_prev_cens)

Observations:
         Total  Left-censored     Uncensored Right-censored 
          2237              6           2231              0 

Coefficients:
                  Estimate Std. error t value  Pr(> t)    
(Intercept)      -0.066292   0.016999  -3.900 9.63e-05 ***
pers_o            0.037621   0.021789   1.727   0.0842 .  
pers_c            0.051433   0.022403   2.296   0.0217 *  
pers_e            0.034091   0.018704   1.823   0.0683 .  
pers_a            0.036618   0.023654   1.548   0.1216    
pers_n            0.036382   0.022615   1.609   0.1077    
airport_distance -0.129827   0.022412  -5.793 6.93e-09 ***
republican       -0.114034   0.020184  -5.650 1.61e-08 ***
medage           -0.124353   0.018030  -6.897 5.30e-12 ***
male             -0.041640   0.021246  -1.960   0.0500 .  
popdens           0.110682   0.016589   6.672 2.52e-11 ***
manufact         -0.013706   0.019420  -0.706   0.4803    
tourism           0.007299   0.021471   0.340   0.7339    
academics        -0.146510   0.030592  -4.789 1.67e-06 ***
medinc            0.172215   0.026205   6.572 4.97e-11 ***
physician_pc     -0.044685   0.017889  -2.498   0.0125 *  
onset_prev       -0.450417   0.020633 -21.829  < 2e-16 ***
logSigma         -0.290067   0.014981 -19.362  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Newton-Raphson maximisation, 4 iterations
Return code 2: successive function values within tolerance limit
Log-likelihood: -2526.32 on 18 Df

tobit_slope_prev_ctrl %>% confint(level=0.9)

                           5 %         95 %
(Intercept)      -0.0942521972 -0.038331633
pers_o            0.0017801165  0.073460914
pers_c            0.0145836752  0.088282659
pers_e            0.0033263545  0.064856355
pers_a           -0.0022888725  0.075525539
pers_n           -0.0008167478  0.073580311
airport_distance -0.1666914394 -0.092961616
republican       -0.1472344136 -0.080833897
medage           -0.1540092486 -0.094697335
male             -0.0765867008 -0.006693169
popdens           0.0833962602  0.137967837
manufact         -0.0456494616  0.018237042
tourism          -0.0280170988  0.042615471
academics        -0.1968297667 -0.096190872
medinc            0.1291117067  0.215317728
physician_pc     -0.0741096257 -0.015260579
onset_prev       -0.4843554389 -0.416477572
logSigma         -0.3147085015 -0.265424820

Predict prevalence slopes with censored regression

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_slope_prev <- cforest(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n + 
                           airport_distance + republican + medage + male + popdens + 
                           manufact + tourism + academics + medinc + physician_pc +
                            onset_prev,
                           data = df_us_slope_prev, 
                         controls = ctrls)

crf_slope_prev_varimp <- varimp(crf_slope_prev, nperm = 1)
crf_slope_prev_varimp_cond <- varimp(crf_slope_prev, conditional = T, nperm = 1)

crf_slope_prev_varimp

          pers_o           pers_c           pers_e           pers_a 
     0.009038452      0.005962953      0.002011857      0.006557738 
          pers_n airport_distance       republican           medage 
     0.005014302      0.015549759      0.017518869      0.028666088 
            male          popdens         manufact          tourism 
     0.007709733      0.146701401      0.015152517      0.003916224 
       academics           medinc     physician_pc       onset_prev 
     0.010309543      0.012499734      0.002298431      0.355694241 

crf_slope_prev_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_slope_prev_varimp_cond

          pers_o           pers_c           pers_e           pers_a 
     0.014171474      0.005819250      0.002895896      0.006566015 
          pers_n airport_distance       republican           medage 
     0.006903428      0.017787026      0.019239684      0.031186235 
            male          popdens         manufact          tourism 
     0.007856779      0.158747164      0.014515458      0.005271982 
       academics           medinc     physician_pc       onset_prev 
     0.014177035      0.014597721      0.003466332      0.399068410 

crf_slope_prev_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

CRF predicting slopes

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_slope_prev <- cforest(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n + 
                           airport_distance + republican + medage + male + popdens + 
                           manufact + tourism + academics + medinc + physician_pc +
                            onset_prev,
                           data = df_us_slope_prev, 
                         controls = ctrls)

crf_slope_prev_varimp <- varimp(crf_slope_prev, nperm = 1)
crf_slope_prev_varimp_cond <- varimp(crf_slope_prev, conditional = T, nperm = 1)

crf_slope_prev_varimp
crf_slope_prev_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_slope_prev_varimp_cond
crf_slope_prev_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

Predict Social Distancing

Change point analysis

# keep only counties with full data
fips_complete <- df_us_socdist_scaled %>% 
  group_by(county_fips) %>% 
  summarize(n = n()) %>% 
  filter(n==max(.$n)) %>% 
  .$county_fips

# run changepoint analysis
df_us_socdist_cpt_results <- df_us_socdist_scaled %>% 
  select(county_fips, socdist_single_tile) %>%
  filter(county_fips %in% fips_complete) %>% 
  split(.$county_fips) %>%
  map(~ cpt.meanvar(as.vector(.$socdist_single_tile),
                    #penalty = 'Asymptotic',
                    class=TRUE,
                    param.estimates=TRUE,
                    Q=1,
                    test.stat = 'Normal'))

df_us_socdist_cpt_results_2 <- df_us_socdist_scaled %>% 
  select(county_fips, socdist_tiles) %>%
  filter(county_fips %in% fips_complete) %>% 
  split(.$county_fips) %>%
  map(~ cpt.meanvar(as.vector(.$socdist_tiles),
                    #penalty = 'Asymptotic',
                    class=TRUE,
                    param.estimates=TRUE,
                    Q=1,
                    test.stat = 'Normal'))

# calculate change point
df_us_socdist_cpt_day <- df_us_socdist_cpt_results %>% 
  map(cpts) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(cpt_day_socdist = '.') %>%
  rownames_to_column('county_fips')

df_us_socdist_cpt_day_2 <- df_us_socdist_cpt_results_2 %>% 
  map(cpts) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(cpt_day_socdist_2 = '.') %>%
  rownames_to_column('county_fips')

# calculate mean differences
df_us_socdist_cpt_mean_diff <- df_us_socdist_cpt_results %>% 
  map(param.est) %>% 
  map(~ .$mean) %>% 
  map(~ .[2]) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(mean_diff_socdist = '.') %>%
  rownames_to_column('county_fips')

df_us_socdist_cpt_mean_diff_2 <- df_us_socdist_cpt_results_2 %>% 
  map(param.est) %>% 
  map(~ .$mean) %>% 
  map(~ -.[2]) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(mean_diff_socdist_2 = '.') %>%
  rownames_to_column('county_fips')

# calculate means
df_us_socdist_mean <- df_us_socdist_scaled %>%
  group_by(county_fips) %>%
  summarize(mean_socdist = mean(socdist_single_tile))

df_us_socdist_mean_2 <- df_us_socdist_scaled %>%
  group_by(county_fips) %>%
  summarize(mean_socdist_2 = -mean(socdist_tiles))

# merge with county data
df_us_cpt_socdist <- df_us_socdist_scaled %>% 
  select(-time, -socdist_single_tile, -socdist_tiles) %>%
  distinct() %>% 
  left_join(df_us_socdist_cpt_day, by='county_fips') %>%
  left_join(df_us_socdist_cpt_day_2, by='county_fips') %>%
  left_join(df_us_socdist_cpt_mean_diff, by='county_fips') %>%
  left_join(df_us_socdist_cpt_mean_diff_2, by='county_fips') %>%
  left_join(df_us_socdist_mean, by='county_fips') %>%
  left_join(df_us_socdist_mean_2, by='county_fips') %>%
  left_join(select(df_us_onset_prev, county_fips, onset_prev), by='county_fips') %>%
  left_join(select(df_us_slope_prev, county_fips, slope_prev), by='county_fips') 

# standardize mean/var differences
df_us_cpt_socdist <- df_us_cpt_socdist %>% 
  mutate(mean_diff_socdist = scale(mean_diff_socdist),
         mean_diff_socdist_2 = scale(mean_diff_socdist_2),
         mean_socdist = scale(mean_socdist),
         mean_socdist_2 = scale(mean_socdist_2))

# handle censored data
df_us_cpt_socdist <- df_us_cpt_socdist %>% 
  mutate(cpt_day_socdist = ifelse(is.na(cpt_day_socdist), as.numeric(diff(range(df_us_socdist$date))), cpt_day_socdist)) %>% 
  mutate(event = ifelse(cpt_day_socdist >= as.numeric(diff(range(df_us_socdist$date))), 0, 1))

df_us_cpt_socdist$cpt_day_socdist %>% hist()
df_us_cpt_socdist$mean_diff_socdist %>% hist()
df_us_cpt_socdist$mean_socdist %>% hist()


df_us_cpt_socdist$cpt_day_socdist_2 %>% hist()
df_us_cpt_socdist$mean_diff_socdist_2 %>% hist()
df_us_cpt_socdist$mean_socdist_2 %>% hist()

# predict hazard from personality
cox_cpt_socdist <- coxph(Surv(cpt_day_socdist, event) ~ 
                           pers_o + pers_c + pers_e + pers_a + pers_n, 
                  data = df_us_cpt_socdist)
cox_cpt_socdist %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n, data = df_us_cpt_socdist)

  n= 2380, number of events= 2379 

            coef exp(coef)  se(coef)      z Pr(>|z|)    
pers_o -0.189086  0.827715  0.023503 -8.045 8.61e-16 ***
pers_c -0.033629  0.966931  0.027365 -1.229  0.21912    
pers_e  0.004178  1.004187  0.023905  0.175  0.86126    
pers_a -0.084190  0.919256  0.028674 -2.936  0.00332 ** 
pers_n  0.060781  1.062667  0.027280  2.228  0.02588 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

       exp(coef) exp(-coef) lower .95 upper .95
pers_o    0.8277     1.2081    0.7905    0.8667
pers_c    0.9669     1.0342    0.9164    1.0202
pers_e    1.0042     0.9958    0.9582    1.0524
pers_a    0.9193     1.0878    0.8690    0.9724
pers_n    1.0627     0.9410    1.0073    1.1210

Concordance= 0.593  (se = 0.007 )
Likelihood ratio test= 121.9  on 5 df,   p=<2e-16
Wald test            = 124.3  on 5 df,   p=<2e-16
Score (logrank) test = 123.9  on 5 df,   p=<2e-16

# predict hazard from personality with controls
cox_cpt_socdist_ctrl <- coxph(Surv(cpt_day_socdist, event) ~ 
                                 pers_o + pers_c + pers_e + pers_a + pers_n + 
                           airport_distance + republican + medage + male + popdens + 
                           manufact + tourism + academics + medinc + physician_pc +
                           onset_prev + slope_prev ,
                  data = df_us_cpt_socdist)
cox_cpt_socdist_ctrl %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n + airport_distance + republican + 
    medage + male + popdens + manufact + tourism + academics + 
    medinc + physician_pc + onset_prev + slope_prev, data = df_us_cpt_socdist)

  n= 2196, number of events= 2195 
   (184 observations deleted due to missingness)

                      coef exp(coef)  se(coef)      z Pr(>|z|)    
pers_o           -0.079525  0.923555  0.029478 -2.698 0.006981 ** 
pers_c           -0.036833  0.963837  0.028968 -1.272 0.203550    
pers_e            0.010726  1.010783  0.025616  0.419 0.675430    
pers_a           -0.086909  0.916761  0.031916 -2.723 0.006468 ** 
pers_n            0.024748  1.025056  0.031361  0.789 0.430049    
airport_distance  0.004603  1.004614  0.026869  0.171 0.863972    
republican        0.010261  1.010314  0.026991  0.380 0.703818    
medage            0.007869  1.007900  0.022040  0.357 0.721057    
male              0.002735  1.002739  0.023670  0.116 0.908011    
popdens          -0.225041  0.798483  0.068056 -3.307 0.000944 ***
manufact          0.099373  1.104478  0.024223  4.102 4.09e-05 ***
tourism          -0.075111  0.927641  0.026289 -2.857 0.004275 ** 
academics         0.123576  1.131536  0.042617  2.900 0.003735 ** 
medinc           -0.033975  0.966596  0.036184 -0.939 0.347760    
physician_pc     -0.015227  0.984888  0.025260 -0.603 0.546629    
onset_prev       -0.004160  0.995849  0.003406 -1.221 0.222018    
slope_prev       -0.241254  0.785642  0.028769 -8.386  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                 exp(coef) exp(-coef) lower .95 upper .95
pers_o              0.9236     1.0828    0.8717    0.9785
pers_c              0.9638     1.0375    0.9106    1.0201
pers_e              1.0108     0.9893    0.9613    1.0628
pers_a              0.9168     1.0908    0.8612    0.9759
pers_n              1.0251     0.9756    0.9639    1.0900
airport_distance    1.0046     0.9954    0.9531    1.0589
republican          1.0103     0.9898    0.9583    1.0652
medage              1.0079     0.9922    0.9653    1.0524
male                1.0027     0.9973    0.9573    1.0504
popdens             0.7985     1.2524    0.6988    0.9124
manufact            1.1045     0.9054    1.0533    1.1582
tourism             0.9276     1.0780    0.8811    0.9767
academics           1.1315     0.8838    1.0409    1.2301
medinc              0.9666     1.0346    0.9004    1.0376
physician_pc        0.9849     1.0153    0.9373    1.0349
onset_prev          0.9958     1.0042    0.9892    1.0025
slope_prev          0.7856     1.2728    0.7426    0.8312

Concordance= 0.627  (se = 0.008 )
Likelihood ratio test= 254.1  on 17 df,   p=<2e-16
Wald test            = 239.4  on 17 df,   p=<2e-16
Score (logrank) test = 244.1  on 17 df,   p=<2e-16

for(i in head(df_us_socdist_cpt_results_2, 5)){
  plot(i)
}

Predicting change points with time-to-event regression

lm_meandiff_socdist <- lm(mean_diff_socdist ~ 
                            pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_us_cpt_socdist)
lm_meandiff_socdist %>% summary()

Call:
lm(formula = mean_diff_socdist ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n, data = df_us_cpt_socdist)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2026 -0.5902 -0.0725  0.5561  5.2119 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -8.612e-16  1.865e-02   0.000 1.000000    
pers_o       3.542e-01  2.022e-02  17.515  < 2e-16 ***
pers_c      -1.084e-01  2.530e-02  -4.284 1.91e-05 ***
pers_e       8.717e-02  2.080e-02   4.191 2.88e-05 ***
pers_a      -9.711e-02  2.583e-02  -3.760 0.000174 ***
pers_n       9.444e-02  2.360e-02   4.002 6.48e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9098 on 2374 degrees of freedom
Multiple R-squared:  0.174, Adjusted R-squared:  0.1723 
F-statistic:   100 on 5 and 2374 DF,  p-value: < 2.2e-16

lm_meandiff_socdist %>% confint(level=0.9)

                    5 %        95 %
(Intercept) -0.03068621  0.03068621
pers_o       0.32088919  0.38743637
pers_c      -0.15000916 -0.06674567
pers_e       0.05294230  0.12139944
pers_a      -0.13960649 -0.05460713
pers_n       0.05560451  0.13326663

lm_meandiff_socdist_ctrl <- lm(mean_diff_socdist ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n + 
                              airport_distance + republican + medage + male + popdens + 
                              manufact + tourism + academics + medinc + physician_pc +
                              onset_prev + slope_prev ,
                            data = df_us_cpt_socdist)
lm_meandiff_socdist_ctrl %>% summary()

Call:
lm(formula = mean_diff_socdist ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + airport_distance + republican + medage + male + 
    popdens + manufact + tourism + academics + medinc + physician_pc + 
    onset_prev + slope_prev, data = df_us_cpt_socdist)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.5265 -0.4524 -0.0495  0.4317  3.7339 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       0.384566   0.157347   2.444 0.014602 *  
pers_o            0.050637   0.019835   2.553 0.010749 *  
pers_c           -0.038501   0.020805  -1.851 0.064371 .  
pers_e            0.030966   0.017349   1.785 0.074418 .  
pers_a           -0.067550   0.022233  -3.038 0.002407 ** 
pers_n            0.193759   0.021160   9.157  < 2e-16 ***
airport_distance -0.053213   0.016944  -3.140 0.001709 ** 
republican       -0.206865   0.018365 -11.264  < 2e-16 ***
medage            0.316357   0.016330  19.373  < 2e-16 ***
male              0.061390   0.016832   3.647 0.000271 ***
popdens           0.329284   0.038020   8.661  < 2e-16 ***
manufact         -0.075423   0.017932  -4.206 2.70e-05 ***
tourism           0.058192   0.018085   3.218 0.001311 ** 
academics         0.161068   0.029716   5.420 6.61e-08 ***
medinc            0.129013   0.025241   5.111 3.48e-07 ***
physician_pc      0.001339   0.016398   0.082 0.934927    
onset_prev       -0.005433   0.002330  -2.332 0.019797 *  
slope_prev        0.151238   0.019733   7.664 2.69e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6911 on 2178 degrees of freedom
  (184 observations deleted due to missingness)
Multiple R-squared:  0.498, Adjusted R-squared:  0.4941 
F-statistic: 127.1 on 17 and 2178 DF,  p-value: < 2.2e-16

lm_meandiff_socdist_ctrl %>% confint(level=0.9)

                          5 %         95 %
(Intercept)       0.125642215  0.643489450
pers_o            0.017998062  0.083276387
pers_c           -0.072737018 -0.004264821
pers_e            0.002417378  0.059515385
pers_a           -0.104134763 -0.030964564
pers_n            0.158939161  0.228579562
airport_distance -0.081095530 -0.025330047
republican       -0.237085304 -0.176644976
medage            0.289485938  0.343228106
male              0.033691376  0.089088738
popdens           0.266719876  0.391847728
manufact         -0.104931822 -0.045915121
tourism           0.028432430  0.087951016
academics         0.112168212  0.209967646
medinc            0.087478364  0.170548176
physician_pc     -0.025644312  0.028322289
onset_prev       -0.009266471 -0.001599000
slope_prev        0.118766109  0.183708968

Linear models predicting mean differences

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_meandiff_socdist <- cforest(mean_diff_socdist ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n +
                                  airport_distance + republican + medage + male + popdens + 
                                  manufact + tourism + academics + medinc + physician_pc,
                           data = df_us_cpt_socdist,
                         controls = ctrls)

crf_meandiff_socdist_varimp <- varimp(crf_meandiff_socdist, nperm = 1)
crf_meandiff_socdist_varimp_cond <- varimp(crf_meandiff_socdist, conditional = T, nperm = 1)

crf_meandiff_socdist_varimp

          pers_o           pers_c           pers_e           pers_a 
     0.023820991      0.025051356      0.001157161      0.057851204 
          pers_n airport_distance       republican           medage 
     0.031209835      0.011115324      0.147250970      0.177245361 
            male          popdens         manufact          tourism 
     0.009491305      0.231710899      0.027431972      0.027842457 
       academics           medinc     physician_pc 
     0.072246678      0.066667134      0.001334095 

crf_meandiff_socdist_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

crf_meandiff_socdist_varimp_cond

          pers_o           pers_c           pers_e           pers_a 
     0.025376206      0.027343070      0.001738617      0.062022318 
          pers_n airport_distance       republican           medage 
     0.033256945      0.010758335      0.154836773      0.182728962 
            male          popdens         manufact          tourism 
     0.009199907      0.238997444      0.034532457      0.030523903 
       academics           medinc     physician_pc 
     0.078720326      0.065719158      0.002819058 

crf_meandiff_socdist_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

CRF predicting mean difference

# predict hazard from personality
cox_cpt_socdist_2 <- coxph(Surv(cpt_day_socdist_2, event) ~ 
                           pers_o + pers_c + pers_e + pers_a + pers_n, 
                  data = df_us_cpt_socdist)
cox_cpt_socdist_2 %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist_2, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n, data = df_us_cpt_socdist)

  n= 2380, number of events= 2379 

           coef exp(coef) se(coef)      z Pr(>|z|)    
pers_o -0.14730   0.86303  0.02283 -6.453 1.10e-10 ***
pers_c -0.11744   0.88919  0.02688 -4.369 1.25e-05 ***
pers_e  0.03018   1.03064  0.02355  1.282  0.19998    
pers_a  0.03223   1.03276  0.02720  1.185  0.23607    
pers_n  0.07865   1.08183  0.02638  2.982  0.00287 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

       exp(coef) exp(-coef) lower .95 upper .95
pers_o    0.8630     1.1587    0.8253    0.9025
pers_c    0.8892     1.1246    0.8436    0.9373
pers_e    1.0306     0.9703    0.9841    1.0793
pers_a    1.0328     0.9683    0.9791    1.0893
pers_n    1.0818     0.9244    1.0273    1.1392

Concordance= 0.598  (se = 0.008 )
Likelihood ratio test= 109.1  on 5 df,   p=<2e-16
Wald test            = 112.2  on 5 df,   p=<2e-16
Score (logrank) test = 111.5  on 5 df,   p=<2e-16

# predict hazard from personality with controls
cox_cpt_socdist_ctrl_2 <- coxph(Surv(cpt_day_socdist_2, event) ~ 
                                 pers_o + pers_c + pers_e + pers_a + pers_n + 
                           airport_distance + republican + medage + male + popdens + 
                           manufact + tourism + academics + medinc + physician_pc +
                           onset_prev + slope_prev ,
                  data = df_us_cpt_socdist)
cox_cpt_socdist_ctrl_2 %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist_2, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n + airport_distance + republican + 
    medage + male + popdens + manufact + tourism + academics + 
    medinc + physician_pc + onset_prev + slope_prev, data = df_us_cpt_socdist)

  n= 2196, number of events= 2195 
   (184 observations deleted due to missingness)

                       coef  exp(coef)   se(coef)      z Pr(>|z|)    
pers_o           -0.0977782  0.9068500  0.0295918 -3.304 0.000952 ***
pers_c           -0.1102061  0.8956495  0.0290631 -3.792 0.000149 ***
pers_e            0.0423876  1.0432987  0.0258637  1.639 0.101238    
pers_a            0.0196323  1.0198262  0.0307933  0.638 0.523766    
pers_n            0.0577244  1.0594230  0.0312407  1.848 0.064641 .  
airport_distance -0.0204996  0.9797090  0.0257576 -0.796 0.426109    
republican       -0.0204758  0.9797324  0.0268872 -0.762 0.446332    
medage            0.0001654  1.0001654  0.0227102  0.007 0.994188    
male             -0.0216697  0.9785634  0.0235641 -0.920 0.357780    
popdens          -0.1960015  0.8220110  0.0638107 -3.072 0.002129 ** 
manufact          0.1517729  1.1638959  0.0251996  6.023 1.71e-09 ***
tourism          -0.0607958  0.9410154  0.0263591 -2.306 0.021086 *  
academics         0.1096653  1.1159045  0.0410448  2.672 0.007544 ** 
medinc            0.0207527  1.0209695  0.0343697  0.604 0.545972    
physician_pc     -0.0434169  0.9575121  0.0247516 -1.754 0.079413 .  
onset_prev       -0.0054441  0.9945707  0.0034311 -1.587 0.112583    
slope_prev       -0.2044823  0.8150692  0.0287518 -7.112 1.14e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                 exp(coef) exp(-coef) lower .95 upper .95
pers_o              0.9068     1.1027    0.8557    0.9610
pers_c              0.8956     1.1165    0.8461    0.9481
pers_e              1.0433     0.9585    0.9917    1.0975
pers_a              1.0198     0.9806    0.9601    1.0833
pers_n              1.0594     0.9439    0.9965    1.1263
airport_distance    0.9797     1.0207    0.9315    1.0304
republican          0.9797     1.0207    0.9294    1.0327
medage              1.0002     0.9998    0.9566    1.0457
male                0.9786     1.0219    0.9344    1.0248
popdens             0.8220     1.2165    0.7254    0.9315
manufact            1.1639     0.8592    1.1078    1.2228
tourism             0.9410     1.0627    0.8936    0.9909
academics           1.1159     0.8961    1.0297    1.2094
medinc              1.0210     0.9795    0.9545    1.0921
physician_pc        0.9575     1.0444    0.9122    1.0051
onset_prev          0.9946     1.0055    0.9879    1.0013
slope_prev          0.8151     1.2269    0.7704    0.8623

Concordance= 0.637  (se = 0.008 )
Likelihood ratio test= 232.8  on 17 df,   p=<2e-16
Wald test            = 227.3  on 17 df,   p=<2e-16
Score (logrank) test = 230.4  on 17 df,   p=<2e-16

Predicting change points with time-to-event regression

lm_meandiff_socdist_2 <- lm(mean_diff_socdist_2 ~ 
                            pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_us_cpt_socdist)
lm_meandiff_socdist_2 %>% summary()

Call:
lm(formula = mean_diff_socdist_2 ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n, data = df_us_cpt_socdist)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.2029 -0.4306  0.0314  0.4752  3.5934 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -5.361e-16  1.890e-02   0.000  1.00000    
pers_o       3.433e-01  2.049e-02  16.753  < 2e-16 ***
pers_c      -1.273e-01  2.564e-02  -4.962 7.46e-07 ***
pers_e       8.647e-02  2.108e-02   4.101 4.25e-05 ***
pers_a      -1.013e-02  2.618e-02  -0.387  0.69876    
pers_n      -6.321e-02  2.392e-02  -2.643  0.00828 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9221 on 2374 degrees of freedom
Multiple R-squared:  0.1515,    Adjusted R-squared:  0.1497 
F-statistic: 84.79 on 5 and 2374 DF,  p-value: < 2.2e-16

lm_meandiff_socdist_2 %>% confint(level=0.9)

                    5 %        95 %
(Intercept) -0.03110172  0.03110172
pers_o       0.30962103  0.37706929
pers_c      -0.16944669 -0.08505576
pers_e       0.05177379  0.12115788
pers_a      -0.05320705  0.03294326
pers_n      -0.10256368 -0.02384996

lm_meandiff_socdist_ctrl_2 <- lm(mean_diff_socdist_2 ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n + 
                              airport_distance + republican + medage + male + popdens + 
                              manufact + tourism + academics + medinc + physician_pc +
                              onset_prev + slope_prev ,
                            data = df_us_cpt_socdist)
lm_meandiff_socdist_ctrl_2 %>% summary()

Call:
lm(formula = mean_diff_socdist_2 ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n + airport_distance + republican + medage + 
    male + popdens + manufact + tourism + academics + medinc + 
    physician_pc + onset_prev + slope_prev, data = df_us_cpt_socdist)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.5540 -0.2750  0.0349  0.3743  2.4846 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       0.797330   0.172285   4.628 3.91e-06 ***
pers_o            0.098013   0.021718   4.513 6.73e-06 ***
pers_c           -0.065237   0.022780  -2.864  0.00423 ** 
pers_e            0.028005   0.018996   1.474  0.14056    
pers_a            0.039602   0.024343   1.627  0.10392    
pers_n            0.094196   0.023169   4.066 4.96e-05 ***
airport_distance -0.080890   0.018553  -4.360 1.36e-05 ***
republican       -0.162802   0.020108  -8.096 9.32e-16 ***
medage            0.101116   0.017880   5.655 1.76e-08 ***
male              0.018963   0.018430   1.029  0.30365    
popdens           0.216815   0.041629   5.208 2.09e-07 ***
manufact         -0.053730   0.019635  -2.736  0.00626 ** 
tourism           0.054977   0.019801   2.776  0.00554 ** 
academics         0.010810   0.032537   0.332  0.73975    
medinc            0.262251   0.027637   9.489  < 2e-16 ***
physician_pc     -0.018855   0.017954  -1.050  0.29375    
onset_prev       -0.011561   0.002551  -4.532 6.16e-06 ***
slope_prev        0.175028   0.021606   8.101 8.99e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7567 on 2178 degrees of freedom
  (184 observations deleted due to missingness)
Multiple R-squared:  0.4256,    Adjusted R-squared:  0.4212 
F-statistic: 94.94 on 17 and 2178 DF,  p-value: < 2.2e-16

lm_meandiff_socdist_ctrl_2 %>% confint(level=0.9)

                           5 %         95 %
(Intercept)       0.5138262725  1.080833598
pers_o            0.0622754890  0.133750793
pers_c           -0.1027232215 -0.027750846
pers_e           -0.0032540685  0.059264346
pers_a           -0.0004560502  0.079660316
pers_n            0.0560700586  0.132321538
airport_distance -0.1114197918 -0.050360399
republican       -0.1958908439 -0.129712816
medage            0.0716937998  0.130537801
male             -0.0113653754  0.049290949
popdens           0.1483115400  0.285317985
manufact         -0.0860393375 -0.021420084
tourism           0.0223922310  0.087561014
academics        -0.0427319068  0.064351788
medinc            0.2167730070  0.307728772
physician_pc     -0.0484002307  0.010689509
onset_prev       -0.0157587671 -0.007363411
slope_prev        0.1394735086  0.210581500

Linear models predicting mean differences

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_meandiff_socdist_2 <- cforest(mean_diff_socdist_2 ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n +
                                  airport_distance + republican + medage + male + popdens + 
                                  manufact + tourism + academics + medinc + physician_pc,
                           data = df_us_cpt_socdist,
                         controls = ctrls)

crf_meandiff_socdist_varimp <- varimp(crf_meandiff_socdist_2, nperm = 1)
crf_meandiff_socdist_varimp_cond <- varimp(crf_meandiff_socdist_2, conditional = T, nperm = 1)

crf_meandiff_socdist_varimp

          pers_o           pers_c           pers_e           pers_a 
     0.037116142      0.011131620      0.001751406      0.015402157 
          pers_n airport_distance       republican           medage 
     0.003470990      0.014201374      0.114029674      0.015711834 
            male          popdens         manufact          tourism 
     0.005168556      0.215045051      0.018685589      0.011515721 
       academics           medinc     physician_pc 
     0.058227126      0.100319631      0.002844632 

crf_meandiff_socdist_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

crf_meandiff_socdist_varimp_cond

          pers_o           pers_c           pers_e           pers_a 
     0.042193802      0.013462952      0.001172980      0.016247157 
          pers_n airport_distance       republican           medage 
     0.004249825      0.015276098      0.116648436      0.015332295 
            male          popdens         manufact          tourism 
     0.005754267      0.220635909      0.019523230      0.014738936 
       academics           medinc     physician_pc 
     0.074979132      0.107775691      0.003803785 

crf_meandiff_socdist_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

CRF predicting mean difference

us_list_results <- list(cox_onset_prev, cox_onset_prev_ctrl, 
     lm_slope_prev, lm_slope_prev_ctrl, 
     cox_cpt_socdist, cox_cpt_socdist_ctrl,
     lm_meandiff_socdist, lm_meandiff_socdist_ctrl,
     cox_cpt_socdist_2, cox_cpt_socdist_ctrl_2,
     lm_meandiff_socdist_2, lm_meandiff_socdist_ctrl_2)

results_names <- list('cox_onset_prev', 'cox_onset_prev_ctrl', 
     'lm_slope_prev', 'lm_slope_prev_ctrl', 
     'cox_cpt_socdist', 'cox_cpt_socdist_ctrl', 
     'lm_meandiff_socdist', 'lm_meandiff_socdist_ctrl',
     'cox_cpt_socdist_2', 'cox_cpt_socdist_ctrl_2',
     'lm_meandiff_socdist_2', 'lm_meandiff_socdist_ctrl_2')

names(us_list_results) <- results_names

save(us_list_results, file="us_list_results.RData")

Export data

us_list_results <- list(cox_onset_prev, cox_onset_prev_ctrl, 
     lm_slope_prev, lm_slope_prev_ctrl, 
     cox_cpt_socdist, cox_cpt_socdist_ctrl,
     lm_meandiff_socdist, lm_meandiff_socdist_ctrl,
     cox_cpt_socdist_2, cox_cpt_socdist_ctrl_2,
     lm_meandiff_socdist_2, lm_meandiff_socdist_ctrl_2)

results_names <- list('cox_onset_prev', 'cox_onset_prev_ctrl', 
     'lm_slope_prev', 'lm_slope_prev_ctrl', 
     'cox_cpt_socdist', 'cox_cpt_socdist_ctrl', 
     'lm_meandiff_socdist', 'lm_meandiff_socdist_ctrl',
     'cox_cpt_socdist_2', 'cox_cpt_socdist_ctrl_2',
     'lm_meandiff_socdist_2', 'lm_meandiff_socdist_ctrl_2')

names(us_list_results) <- results_names

save(us_list_results, file="us_list_results.RData")

write_csv(df_us_slope_prev, '/Users/hp2500/Google Drive/STUDY/Columbia/Research/Corona/Delivery/df_us_slope_prev.csv')
write_csv(df_us_cpt_socdist, '/Users/hp2500/Google Drive/STUDY/Columbia/Research/Corona/Delivery/df_us_cpt_socdist.csv')